function [x] = matlab_fast_solve()
clc; clear;

A = [];
AA = [];
B = [];
bb = [];
D = [];
gdim = [];
maxit = [];
PI = [];
restart = [];
S = [];
tgdof = [];
tol = [];
vgdof = [];
elem = [];




load('/Users/chen/git/mixed-hzelem/test/data/fast_data.mat')

elem = elem+1;

SL = tril(S);
SU = triu(S,1);
SLT = SU';
SUT = SL';


amg_option = struct('solver','NO','solvermaxit',3,'preconditioner','V','printlevel',0,'setupflag',true,'N0',500);


b = ones(size(A,1),1);
%[~,~,Ai,Bi,BBi,Res,Pro,cl] = amg(A,b,amg_option);
%[~,~,Ai,Bi,BBi,Res,Pro,cl] = mg(A,b,elem,amg_option);
[~,~,Ai,Bi,BTi,Res,Pro,isFreeDof] = mg(A,b,elem,amg_option);

amg_option.solver = 'CG';
amg_option.setupflag = false;


D = D';



tic;
[x,~,relres,iter,resvec] = gmres(AA,bb',restart,tol,maxit,@linear_operator);
fprintf('fast solve time: %f\n',toc)
fprintf('inter number: %f\n',iter)
fprintf('relres: %.3e\n',relres)
%plot(log10(resvec/norm(bb,2)))


    function r = linear_operator(r)
        r(1:tgdof) = r(1:tgdof)./D;
        r(tgdof+1:end) =  r(tgdof+1:end) - B*r(1:tgdof);
        u1 = zeros(vgdof,1);
        for i = 1:10
            u1 = SL\(r(tgdof+1:end)-SU*u1);
        end

        r3 = r(tgdof+1:end)-(SL*u1+SU*u1);

        svgdof = vgdof/gdim;
        for i = 1:gdim
            %u1((i-1)*svgdof+1:i*svgdof) = u1((i-1)*svgdof+1:i*svgdof)+PI*(A\(PI'*r3((i-1)*svgdof+1:i*svgdof)));
            %u1((i-1)*svgdof+1:i*svgdof) = u1((i-1)*svgdof+1:i*svgdof)+PI*amg(A,PI'*r3((i-1)*svgdof+1:i*svgdof),amg_option,Ai,Bi,BBi,Res,Pro,cl);
            u1((i-1)*svgdof+1:i*svgdof) = u1((i-1)*svgdof+1:i*svgdof)+PI*mg(A,PI'*r3((i-1)*svgdof+1:i*svgdof),elem,amg_option,Ai,Bi,BTi,Res,Pro,isFreeDof);
        end

        for i = 1:10
            u1 = SUT\(r(tgdof+1:end)-SLT*u1);
        end


        r(1:tgdof) = r(1:tgdof)+(B'*u1)./D;
        r(tgdof+1:end) = -u1;
    end



    function u = linear_operators(r)
        u = zeros(tgdof+vgdof,1);
        u(1:tgdof) = r(1:tgdof)./D;
        u(tgdof+1:end) = S\r(tgdof+1:end);
        %u(tgdof+1:end) = pcg(S,r(tgdof+1:end),1e-6,maxit,L1,L1');
   

        
    end



end



